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An extremely easy method for accurately calculating the 
adiabatic connection of density functional theory is presented, 
and its accuracy tested on both Hooke's atom and the He 
atom. The method is easy because calculations are needed 
only for different values of parameters in the external poten- 
tial, which can be achieved with almost any electronic struc- 
ture code. 



I. INTRODUCTION 

Density functional theory has become a popular com- 
putational method in quantum chemistry, because of its 
ability to handle large molecules accurately but relatively 
inexpensively ||l|,|[. This success is based on the avail- 
ability of reliable accurate approximate functionals, and 
there is a constant need for still further accuracy. The 
goal of atomization energy errors being reliably less than 
1 kcal/mol has not yet been achieved. 

An important step forward in this search for accuracy 
came when Becke mixed a fraction of exact exchange 
with a generalized gradient approximation (GGA), and 
reduced errors by a factor of two or three [p|. Such hy- 
brid functionals, e.g., B3LYP, are now in common use, 
but their underlying justification comes from the adia- 
batic decomposition of density functional theory. Ini- 
tially, the mixing parameters used were determined em- 
pirically. Later, it was shown that these parameters could 
be derived non-empirically [0Jq] , and that a single univer- 
sal mixing coefficient (25%) could be rationalized based 
on the performance of MP theory for molecules 0] . Most 
recently, new functionals have been proposed which use 
this adiabatic decomposition in much detail 0. 

The adiabatic decomposition of an electronic system 
is very simple conceptually. Imagine multiplying the 
electron-electron repulsion by a coupling constant A. 
Now imagine varying A, while keeping the electron den- 
sity p(r) fixed. This differs from traditional perturbative 
methods, e.g., MoUer-Plesset [g[, because the external po- 
tential must be altered at each A to keep the density fixed. 
At A = 1, we have the physical, interacting electronic 
system. But as A is reduced to zero, keeping the density 
fixed, the electrons become those of the non-interacting. 



Kohn-Sham system, and the potential morphs into the 
Kohn-Sham potential. All Kohn-Sham DFT calculations 
are actually performed on this non-interacting system, 
and the physical ground-state energy deduced from it. 
The adiabatic connection provides a continuous connec- 
tion between the interacting system and its Kohn-Sham 
analog. 

The only part of the total energy which must be 
approximated in such a Kohn-Sham calculation is the 
exchange-correlation energy, Exc[p\, as a functional of 
the (spin) density. A further value of the coupling con- 
stant is that, through the Hellmann-Feynman theorem, 
this energy can be written as an intergal over the purely 
potential contribution [pHlCJ]: 



Sxc = / dAf/xc(A), C/xc(A) = (*lKo|*^) - U, (1) 
Jo 

where T4o is the Coulomb interaction between electrons, 
^'^ is the wavefunction at coupling constant A, and U 
is the Hartree electrostatic energy. This integral is the 
adiabatic connection formula. Hybrid functionals are 
based on the fact that this integrand, when applied to 
the exchange-correlation contribution to a bond dissoci- 
ation energy, is often not well-approximated by GGA's, 
due to their lack of static correlation |^,Q. This can 
be partially corrected for most molecules by mixing in a 
fraction of exact exchange. Construction of functionals 
based on this insight is often referred to as the adiabatic 
connection method (ACM). 

Thus accurate approximation of this adiabatic connec- 
tion curve is extremely important to further progress in 
construction of approximate functionals, and benchmark 
cases for small systems are always of interest and help 
in this endeavor. However, this is, in principle, a very 
demanding task. For each value of A, one must solve the 
interacting electronic problem many times in order to 
find the external potential which reproduces the A = 1 
density. Almbladh and Pedroza |jl3] made early heroic 
attempts, but not with the accuracy of modern calcula- 
tions. In the last several years, as the importance of these 
curves has become apparent, several groups in different 
areas have performed adiabatic decomposition calcula- 
tions. Hood et al. jl4| have calculated the curve for bulk 
Si. Colonna and Savi n p5[ have used the slightly differ- 
ent Harris and Jones ]16|| decomposition for both the He 
and Be isoelectronic series. Joubert and Srivastava IT^l 
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FIG. 1. Exact UxcW curve for the k=l/4 Hooke's atom, 
(solid line) and within the PBE correlation functional (dashed 
line), calculated as described in the text (atomic units). 



have used Hylleras-type wavefunctions to calculate these 
curves for the He isoelectronic series. 

All these methods require solving the interacting elec- 
tronic problem with a potential that differs from the orig- 
inal external potential, e.g., —Z/r for atoms. This typ- 
ically makes them difficult to transfer to other systems, 
radically different from the original one, in which other 
codes and approximations are being used. In this pa- 
per, we show how to construct the adiabatic connection 
formula for any atom accurately by doing calculations 
simply for different values of Z, so that no modification 
need be made to an existing wavefunction code. A stan- 
dard calculation is simply run for several different nuclear 
charges, running from the physical value up to Z = cxd. 
The main focus of this paper is to demonstrate the ac- 
curacy of the method. The solid line in Figure 1 is an 
essentially exact curve for C/xc(A) for Hooke's atom (two 
electrons in a harmonic potential), for a spring constant 
of fc = 1/4 (atomic units), calculated from a series of ex- 
act calculations for force constants greater than 1/4, but 
never requiring any calculation with a different external 
potential. The method can be immediately applied to 
accurate calculations for larger atoms, and is currently 
being explored for molecules. 

The paper is divided into several sections. In the next 
section, we outline the basic theory behind our calcula- 
tions. Following that, we present results, both using the 
exact functional, and within a GGA, for both Hooke's 
atom and the He atom. Hooke's atom consists of two 
electrons attached by springs to a center, but interacting 
via a Coulomb repulsion. We summarize our findings in 
the last section. Atomic units (e^ = ?i = me — 1) are 
used throughout, and only spin-unpolarized systems are 
discussed. 



II. THEORY 



A. Formally exact results 



For convenience, we begin with a review of the adia- 
batic connection formalism. We define ^'^[p] to be that 
wavefunction which has density p(r) and minimizes 



Then the kinetic-correlation energy at A is 



(2) 



(3) 



where Ts[p] is the non-interacting Kohn-Sham kinetic 
energy of density p. The potential contribution to 
exchange-correlation at coupling constant A is, 



while the exchange-correlation energy at A is 






r;: 



(4) 



(5) 



For A = 1, the system becomes the true interacting 
quantum-mechanical system, and energies without super- 
scripts refer to A = 1, e.g., E^c = ^xc^- 

All quantities at coupling constants different from one 
are simply related to their full coupling strength counter- 
parts, but evaluated on a scaled density fls]. The most 
important examples are the wavefunction: 



and the exchange-correlation energy 

E^^[p]^X^Exc[pi/xi 



(6) 



(7) 



7^p(7r), 



where P7(r) — 

and \I'-y(ri..rjv) = 7'^/^\I'-y(7ri..7rjv). Thus knowledge 
of how a quantity varies as the density is scaled implies 
knowledge of its coupling constant dependence. 

Quantities evaluated on the Kohn-Sham wavefunction 
vary in a simple way, such as the non-interacting kinetic 
energy Ts: 



T,^[p]=T,[p], or rsK]=7'Ts[p], 
the Hartree electrostatic energy, 

U^[p]^XU[p], or U[p^]^jU[p], 
and the exchange energy 

E^[p]^XEx[pl or Ex[p^]^-fEx[p]. 



(8) 



(9) 



(10) 



Correlation is more sophisticated. Note, however, that 
knowledge of any quantity, Ec[p.y],Tc[P'y], or C/c[P7] as 
a function of 7 for 7 between 1 and 00, i.e., scaling to 



the high density Hmit, is sufficient to determine the adia- 
batic connection for any of them, for A between and 1. 
The most well-known relation is to extract the kinetic- 
correlation piece from Ec [p] |l8[ : 



Tc[p] = -Ea[p] 



dEc [p-, 



d'y 



7=1 



(11) 



To generalize this result to p^, apply Eq. (hW) to p^, and 
make a change of variables in the derivative, to find: 



Tc[p^] = -Ec[p^] +7 



dEc [p-, 
d'y 



(12) 



This can be considered a first-order differential equation 
in 7 for Eclp^y]- Solution of this equation yields 



EciPf] = -7 



^^ V r 1 

■yjTcipy], 



(13) 



where we have used the fact that -Ec[P7] has vanished 
as 7 -^ ||l^. Via Eq. @), these can be turned into 
coupling-constant relations: 



and 



T^ ^Ei. 



Ei = -X 



dE^ 



dX' 
A^ 



T. 



A' 



(14) 



(15) 



The latter is called Bass' relation ^^. Similarly, we can 
extract Uc from Ec[p-y], since Uc — Eq — T^- Thus Eq. 
(|l|) leads to 



Uc[p^]^2Ec[p^] 
which, inverted, is 



dEc[p^ 
d'y 



Ec[Pi 



1 






(16) 



(17) 



Combining Eqs. (Iq) and (O) then leads to simple rela- 
tions between Uc and Tc- Lastly, the coupling-constant 
relation that follows from Eq. (^7|) is 



Ec[p]^j^ ^U'o[pl 



(18) 



and C/xc(A) of Fig. 1 is just the integrand when this 
expression is applied to both exchange and correlation, 
i.e., U^c/X- We emphasize that all these relations follow 
from the well-known Eqs. (0) and (pi]). 



B. Highly accurate approximations 

We begin this section by noting how, when some pa- 
rameter in an external potential is altered, the density 
changes scale, but often does not change shape very 
much. For example, for the two-electron ion, going from 
Z=2 to Z=4 will roughly multiply the density by 8, and 
reduce its length scale by a factor of 2. We use this fact to 
very accurately approximate scaling the density. Then, 
through the relations derived above, we can convert this 
into the coupling-constant dependence. 

Let p(r) be the density of the system we are interested 
in. Suppose we alter some parameter in the external 
potential, and solve the interacting problem, finding some 
density p'(v). If we want to treat p' as an approximation 
to pj, we must first choose a criterion for determining 7. 
In fact, either of Eqs. ^ and (^) could be used, since 
both would be satisfied exactly if p' were truly a scaled 
density. For definiteness, we choose 



A = l/7 = C/[p]/C/[p']. 



(19) 



We can then consider Ec[p'] ~ £^0(^7], Tc[p'] « 7c [^7], 
etc., which we call the bare estimates. As we shall show 
in the results section, these energies usually yield quite 
accurate approximations to the exact quantities. How- 
ever, since we are not working with the true scaled den- 
sity, these energies do not satisfy the relations involving 
scaling derivatives above, such as Eq. (|l^). 

To make a better estimate, we note that, in the case of 
the correlation energy, we can calculate the leading cor- 
rection to Ec[p-y], using the correlation potential. This 
can be constructed by finding the exact Kohn-Sham po- 
tential which generates the density, and subtracting from 
it the external, Hartree, and exchange contributions. 
Then, 



Ec[p,] « Ec[p'] + J d\vc[p']{r) {p,{r) - p'(r)) + 0{Sp)' 

(20) 

As will be shown below, inclusion of this correction leads 
to extremely accurate adiabatic connection curves. We 
generate Tc[p-y] using Eq. (p^), and then Uc = Ec ~ Te- 
la particular, as 7 — > 1, p' — > p, and Eq. (EQ) becomes 
exact to first order in the difference between p' and p. 
Thus Eq. (O) is satisfied exactly. 

Lastly, we apply the same principles to Tg , as a test of 
the closeness of the approximate density, p' . The Euler 
equation for the Kohn-Sham system says that 



5p{r) 



M-Ws(r), 



(21) 



where /i is a constant, and Ws(r) is the Kohn-Sham po- 
tential. Thus, we can also correct the bare Tg estimate, 
to 



Ts[p,] = Up'] ~Jd\vsW]{r) {p,{r) - p'(r)) + 0(6 p)\ 

(22) 



III. RESULTS 
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FIG. 2. Densities used for A = 1 (least p(0)), 0.75, 0.5, 
0.25, and (largest p(0)). 

We present results for two systems, Hooke's atom and 
the He atom. The first of these contains two electrons in 
a harmonic potential, interacting via a Coulomb repul- 
sion. This provides a valid test case for density functional 
theory, because of the Coulomb repulsion between elec- 
trons. It is also an easy system to perform calculations 
on, because the center-of-mass and relative coordinates 
separate, leaving only a one-dimension differential equa- 
tion to solve numerically pi[ . Even this equation has an 
analytic solution for force constant k — 1/4 [2^], which 
is the system on which we will demonstrate most of our 
results. 

To implement our method for the k = 1/4 Hooke's 
atom, we run many calculations at different values of 
k > 1/4, up to fc = 10^, using the method of Ref. [|l|. 
For each calculation, we scale the density so that it looks 
like the original k — 1/4 density, using Eq. ([l9| ) to define 
the scale factor. Several such densities are shown in Fig. 
2, which illustrates how close these densities are. Note 
that the largest error is for A = 0, the non- interacting 
Gaussian density. 

Figure 3 is a plot of the three quantities Ec[pi/x\, 
Tc[pi/\], and Uc[pi/\], all found using the PEE corre- 
lation functional. We choose these quantities to plot, as 
they remain finite in the range A = to 1. For each 
quantity, there are two curves. The dashed curve is the 
bare result of the calculations at different values of the 
force constant, i.e., using E[p'] alone. The solid curve 
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FIG. 3. Correlation energies for the k=l/4 Hooke's atom, 
calculated both with (solid lines) and without (dashed lines) 
the correction term, using the PBE correlation functional. 
The solid lines are indistinguishable from exact results, cal- 
culated by scaling the PBE correlation functional. 



for Ec [pi/7] includes the correction due to the potential, 
in accordance with Eq. (|20|). The f/^ and T^ curves are 
then extracted from this one, using Eqs. ( |l^ ) and (|lq). 
Both the value and first derivative at A = 1 are exact for 
this curve. What is more remarkable is that these curves 
coincide essentially exactly with the correct result, which 
we can deduce by directly applying Eq. (0) to the PBE 
correlation energy functional. Careful numerical calcula- 
tions indicate that the maximum error in our curve is at 
A = 0, and is less than 3xl0~^ Hartrees. This is due to 
the similarity of densities, as shown in Fig. 2, leading to 
very small corrections. 

We also tested other possible prescriptions for choos- 
ing A, such as from the square root of the ratio of non- 
interating kinetic energy densities, as in Eq. (H). This 
gives values for A very close to those of Eq. (||) , and leads 
to no measurable change in our estimate for _Ec[/5i/a]- 

Figure 4 is the analog of Fig. 3, but now calculated 
using exact correlation energies and potentials. Based 
on the remarkable accuracy of Fig. 3, we claim these 
curves are essentially exact. In fact, the change due to 
the potential correction is much smaller for the exact case 
than for PBE, suggesting that the error made in Eq. (|2(]| ) 
will also be smaller. The adiabatic curve of Fig. 1 was 
derived from these curves, since U^ci}^ = ''^t^xc[pi/A]- 

We also ran PBE calculations for lower densities, where 
the shape of the density can change significantly with fc. 
We find that even as low as fc = 10~^, there is a maxi- 
mum error of only 1 millihartree in the corrected Eq [pi/a] 
curve. Beyond this point, the reliability of our method 
might be questioned. However, defining an average r^ 
value via Eq. (6a) of Ref. |g3[, at this point (r^) = 19, 
whereas for fc = 1/4, (r^) is 2.07. Thus for common val- 
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FIG. 4. Correlation energies for the k=l/4 Hooke's atom, 
calculated both with (solid lines) and without (dashed lines) 
the correction term. 



ues of the density, the errors in our procedure are minute. 
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FIG. 5. Non-interacting kinetic energy as a function of cou- 
pling constant. Exact quantity (solid line) is independent of 
A. Dashed line is bare estimate of Ts, the dotted line includes 
the correction. 

As a final test of our method, we return to the non- 
interacting kinetic energy. If the scaHng were exact, the 
non-interacting kinetic energy would scale quadratically, 
as in Eq. (^). In Figure 5, we plot Ts[p^]/7^ exactly, 
approximating Ts[p^] by T's[p'], and including the correc- 
tion of Eq. (p2|). Note that the maximum absolute error, 
at A = 0, is only about -1.6%, in the bare estimate. The 
correction makes the derivative with respect to A exact 
as A — > 1, and overall reduces the error by about a factor 
of 5, to about -0.3%. 

Finally, we discuss our results for the two-electron ion 
series, simply to demonstrate that there is nothing spe- 
cial about Hooke's atom which makes these techniques 
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FIG. 6. Correlation energy for the He atom, calculated 
both without (diamonds) and with (crosses) the correction 
term, using the PBE correlation functional. The solid line 
is the exact result, calculated by scaling the PBE correlation 
functional. 



accurate. We may make bare estimates of the adiabatic 
curves directly from energy data already in the literature. 
Thus by calculating A by Eq. (|9|) using the data in Table 
I of Ref . |2J] , and assigning to each A the bare Uc value 
from the table, we find the corresponding adiabatic con- 
nection curve integrates to Eq = -41.5 mH, as compared 
with an exact value of -42.1 mH for He. We can similarly 
use the Tc data in Eq. (|l^) to get an estimate of -41.2 
mH, so that the difference between these two results is 
a good indicator of the error in both of them. Again, 
we can calculate T^ , and find the largest error at A = 0, 
with a value of 2.916, as compared to the exact value of 
2.867. 

In Fig. 6, we plot PBE calculations of -Ec[pi/a] for the 
He atom, both exactly and using the fake adiabatic con- 
nection formula, both with and without the correction 
term. We use the highly accurate densities of Umrigar 
and Gonze P5[ . Unfortunately, not enough data points 
are available to make smooth plots for this system. Al- 
ready at Z = 3, A = 0.62, so that results for non-integer 
values of Z (between Z — 2 and 3) are needed. How- 
ever, enough data is present to see the clear correction 
the potential makes. We note several interesting features 
of this curve. First, the dependence on scaling is much 
less for He than for the k = 1/4 Hooke's atom. Second, 
the PBE curve contains a minimum, so Eq. (12) implies 
Ec + Tc becomes positive for this scaled density. While 
it has never been rigorously proven, every known case of 
Eq ■\- Tq is negative ||l^,^ . Thus this is probably a (very 
slight) limitation of PBE, which does not occur in PW91, 
and may be related to the difference at Tj = in Fig. 7 
of Ref. [M. Third, the potential correction even picks 
up this feature, and still only has errors of a fraction of 
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FIG. 7. Correlation energy for the He atom, calculated 
both without (diamonds) and with (crosses) the correction 
term. 



7 = 1/A ^ cxD. This effect causes noticeable errors only 
for A < 0.2 in Ec[pi/x]^ but not in Uc or Tc- It also 
means that neither Uc[p\ or T(;:[p], as derived from the 
resulting Ec[p^] curve via Eqs. (Qq) and (p^, is exact at 
A = 1. However, these differences are of the order of 0.1 
millihartree. 

Finally, we are currently investigating if this method 
can be used to calculate accurate adiabatic connection 
curves for molecules [T^ . 
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Figure 7 repeats Fig. 6, but for the exact case. Here 
the corrected curve is much straighter than for FEE, and 
has no minimum. The lines are drawn merely to aid the 
eye. Note that the value of Ec [p] , which is the high- 
density limit of Ec, drops about 1 millihartree, due to 
the difference in shape between the He atom and a pure 
exponential decay (bare estimate). We find Eq = -47.9 
millihartree, in reasonable agreement with the value of 
Colonna and Savin ]13] (-47.5 from Table VI), of Joubert 
and Srivastava p7| (-47.2 from ap in Table II), and of 
Engel and Dreizler ||2| (-48.2 from Table V). The value 
-50.3, reported in Ref. |29], was evaluated on the PBE 
selfconsistent density. 



IV. CONCLUSIONS AND IMPLICATIONS 

We have shown how, with accurate ground-state re- 
sults as a function of the external potential, an accurate 
adiabatic connection curve can be calculated for both 
Hooke's atom at A: = 1/4 and the He atom. We see 
no reason why similar results could not be obtained for 
larger atoms, especially those for which the Kohn-Sham 
potential has already been calculated [pO|. There also 
exist methods for isolating the correlation potential [pl| . 

In the event that the exchange contribution cannot 
be easily isolated, the scheme can still be applied, but 
now to the combined exchange-correlation energy at each 
value of the external parameter. The correction is now 
evaluated using the exchange-correlation potential. On 
Hooke's atom, this yields results almost as accurate as 
those described for the correlation energy alone. Impor- 
tant differences are in the regime A — > 0, as here the 
exchange contribution to the correction blows up, since 
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